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Abstract Territoriality is a phenomenon exhibited throughout nature. On the individ¬ 
ual level, it is the processes by which organisms exclude others of the same species 
from certain parts of space. On the population level, it is the segregation of space into 
separate areas, each used by subsections of the population. Proving mathematically 
that such individual-level processes can cause observed population-level patterns to 
form is necessary for linking these two levels of description in a non-speculative way. 
Previous mathematical analysis has relied upon assuming animals are attracted to a 
central area. This can either be a fixed geographical point, such as a den- or nest-site, 
or a region where they have previously visited. However, recent simulation-based 
studies suggest that this attractive potential is not necessary for territorial pattern for¬ 
mation. Here, we construct a partial differential equation (PDE) model of territorial 
interactions based on the individual-based model (IBM) from those simulation stud¬ 
ies. The resulting PDE does not rely on attraction to spatial locations, but purely on 
conspecific avoidance, mediated via scent-marking. We show analytically that steady- 
state patterns can form, as long as (i) the scent does not decay faster than it takes the 
animal to traverse the terrain, and (ii) the spatial scale over which animals detect scent 
is incorporated into the PDE. As part of the analysis, we develop a general method 
for taking the PDE limit of an IBM that avoids destroying any intrinsic spatial scale 
in the underlying behavioral decisions. 
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1 Introduction 


Territoriality is a wide-spread phenomenon throughout nature. A terr itory is an a rea 
of space used exclusively by an organism, or a group of organisms (IBurll Il943h . It 
is formed by deliberately excluding others of the same species (call ed conspecifics ) 
from the area, either by aggressive confrontations or mutual consent ( Adamsl 2001 ). 
In the last two decades, there have been a number of studies that show analytically 
how territorial patterns can form from the movements and interactions of animals 
( Lewis & Murray . 1993 : Moorcroft & Lewis . 2006t Potts & Lewis . 20141) . These use 
mean-field approximations to model the animals’ behavioral decisions as partial dif¬ 
ferential equations (PDEs), and so enable terri tory formation to be analyzed using 
standard tools from PDE theory ( Murray . 2002h . 


_ Despite their success in uncovering drivers behind space use patterns (iMoorcroft et al. 

20061) . previous analytical models assume an attractive potential influencing the an¬ 
imals’ m ovements. This cou ld either be fidelity to a central place such as a den- or 
nest-site ( Lewis et al . 19971). or a tendency to move towards places that the animal 
has previously visited ( Briscoe et al. . 20021) . However, it is not clear that suc h an 
attractive potential is in fact necessary for territory formation ( Moorcroft 2012 ). 

In this paper, we present a PDE model of territorial pattern formation based purely 
on conspecific avoidance, with no attractive potential. It is based on an individua l 
based model (IBM) of so-called territorial random walkers ( Giuggioli et all 201 la ). 
Previous work used simulation analysis to de monstrate empirically that territories can 
form in this system ( Giuggioli et al. . 201 lal) . Here, we show analytically the circum¬ 
stances under which territorial patterns may form. Specifically, necessary conditions 
for territorial pattern formation include 

- spatial aversion to scent marks 

- scent marks that persist for longer than it takes the animal to traverse the terrain, 

and 

- a reaction to conspecific scent averaged over a small area around the animal. 

As is often the case in ecological applications, it is important that the disc rete 
spatial nature of each interaction is present in the model (IDurrett & LevinLll994l) . In 
the case of territorial interactions, this discreteness is inherent in the fact that animals 
have a non-zero perceptive radius for determining the presence of scent. As part of 
this study, we develop a limiting procedure that enables the transition from IBM to 
PDE without losing this important aspect of spatially discrete interactions. This has 
the pot ential for general use, as previous limiting procedures have often failed in this 
regard ( Durrett & Levin . 1994 ). 

The paper is organized as follows. Section |2] derives the PDE from the IBM 
model. Sections [2 and |4] investigate the conditions under which patterns may form. 
Section|5]gives some concluding remarks. 
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3 ) Agent 1 Agent 2 



Territory of agent 1 Territory of agent 2 


b) Agent 1 Agent 2 



Fig. 1 Pictorial representation of the underlying model. The territory of each agent represents the sites 
containing that agent’s scent. In panel (a), agent 1 is unable to move to the right in the next step, since there 
is active scent of agent 2 there. However, agent 2 can move in either direction. In panel (b), we show the 
case where the lattice spacing, a, is halved, so response to scent is averaged over several sites, given by the 
grey ovals. As the lattice spacing is reduced by a factor of h{a), so the response to scent is averaged over 
2h(a) — 1 sites. In this situation, agent 1 has a higher probability of moving left than right, while agent 2 
has equal probability of moving in both directions. 


2 From the individual-level description to a system of PDFs 

2.1 Description of model 


The individual based model (IBM) is based o n a ID model of territoriality which was 
recently pro posed by Giu 2 doli_ef_fl^ ( 20 11a|i. but then sligh tly modified and studied 
in detail bv iGiuggioli et all (1201 Ibh and IPotts et al.\ (1201211 . The model consists of 
two agents moving on a ID lattice. The agents represent either a single individual 
responsible for territorial defense, or a group of individuals moving together, such as 
a pack or a flock. For example, the former is appropriate when modelling fox {Vulpes 
vulpes) b ehaviour whe re the dominant male in each group marks and secures the 
territory ( Harrisl 1980h . whereas the latter may b e more appropriate for modelling 
wolf (Canis lupus) packs ((Lewis & Murravlll993h . 

Agents move as discrete-time discrete-space neares t-neighbour random wal kers, 
depositing scent marks as they move. In the model of Giuggioli et al. ( 201 Ibh . the 
scent remains present for a finite amount of time, called the active scent time and 
denoted by the symbol TXs ■ Once this time is up, provided the lattice site has not been 
re-scented, the mark is no longer considered by conspeciflcs to be ‘active’. Agents 
cannot move into any lattice site that contains the active scent of another agent (Figure 
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Table 1 Glossary of symbols. The first column shows the symbol, the second a definition, and the third 
whether it pertains to the discrete (lattice) model or the continuous limit or both. Note that some symbols 
are used either as dimensional quantities or their dimensionless equivalents, depending on the context (see 
section 2.5). 


Symbol 

Definition 

Model 

n 

Arbitrary lattice site 

Discrete 

m 

Arbitrary time step 

Discrete 

Ei(n,m) 

For animal i, the probability that there is conspecific scent at (n,m) 

Discrete 

T 

Length of a single time step 

Discrete 

a 

Lattice spacing 

Discrete 

K{a) 

Probability that scent is deposited when the agent visits a lattice site 

Discrete 

h{a) 

Number of lattice sites constituting the agent’s perceptive radius 

Discrete 

frinW) 

Probability of agent i moving to n next jump, given it is at n' at timestep m 

Discrete 

U {n,m) 

Probability that agent 1 is at n at timestep m 

Discrete 

V (n,m) 

Probability that agent 2 is at n at timestep m 

Discrete 

P(n.m) 

Probability that scent of agent 1 is present at n at timestep m 

Discrete 

Q(n,m) 

Probability that scent of agent 2 is present at n at timestep m 

Discrete 

M 

Mean rate of scent decay 

Both 

A 

Mean scent deposition over a unit of space in a unit of time 

Both 

X 

Arbitrary position in continuous space 

Continuous 

t 

Arbitrary continuous time 

Continuous 

u(x.t) 

Probability density function of agent 1 at time t 

Continuous 

v[x.t) 

Probability density function of agent 2 at time t 

Continuous 


Probability that scent of agent 1 is present at (x,t) 

Continuous 

q{x,t) 

Probability that scent of agent 2 is present at (x,?) 

Continuous 

D 

Diffusion constant 

Continuous 

S 

The agent’s perceptive radius 

Continuous 

Pi^d) 

Mean of p(x,t) in a 5-ball around x 

Continuous 

q{x,t) 

Mean of q{xj) in a 5-ball around x 

Continuous 

L 

Width of terrain 

Continuous 

m 

Dimensionless composite variable pLLjX 

Continuous 

e 

Dimensionless composite vaiiable D/LX 

Continuous 


Our model set-up wi ll take three stage s. Stage 1 uses the formalism of coupled 
step selection functions (PoUs_ef_^, 20141) to de scribe a stochastic IBM algorithm 
which generalizes that of Giuggioli et al. ( 201 Ibl) . Stage 2 describes how to derive a 
mean-field probabilistic model from the IBM. Stage 3 involves taking the PDE limit 
of the probabilistic model. 


2.2 Stochastic algorithm for the individual based model 

If unconstrained by scent marks, an agent is simply a nearest-neighbor random walker. 
Therefore the probability that agent i (i G {1,2}) moves from site n' to n is 0,(n|n') = 
l/2if |n —n'l = 1 and 0,(n|n') = 0 otherwise. This function 0,(n|n') is the environment- 
independent movement kernel. 

Now we add the effect of scent marks, which for this paper are thought of as 
constituting of the animal’s ‘environment’. For each agent i, let the environment, 
Ei{n,m), be the probability that there is conspecific scent at lattice site n and timestep 
m. We give two possible definitions for Ei{n,m -f 1), denoted by Ej{npn -(- 1) and 
Ef{n,m-\-l), and both defined in terms of the state of the system at timestep m. The 
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first is given by 


1 if an agent j ^ i is at position n at any time 
£■/ (n, m + 1) = < between m — Tas + 1 and m, 

0 otherwise. 


( 1 ) 


If E} {n,m) — 1 then t here is conspecific scent pre sent, otherw i se the re is not. This is 


the definition used by Giuggioli et al. ( 201 Ibl) and Potts et al. ( 2012 ). 

An alternative to equation ([T]) is the following definition 

{ 1 — jXT with probability K{a), if an agent j ^ i is at 

position n at timestep m, 

{1 — HT)Ef{n,m) otherwise, 

( 2 ) 

where T is the length of a timestep and K{a) is the probability that scent is deposited 
when the animal visits a lattice site. Notice that scent left at timestep m has a proba¬ 
bility 1 — jj.T of remaining present at timestep m+1. 

Introducing K{a) allows us to change the lattice spacing a without changing the 
average distance moved between scent depositions, by insisting that a/K{a) is kept 
constant. From Section 1231 onwards, we will use equation (|2]i to describe scent de¬ 
position and decay. However, the stochastic algorithm of this section can be defined 
equally well using either equation ([T]l or (|2]l. 

We now define the interaction term, which denotes how the scent affects the 
agent’s movement. Animals typically have a perceptive ra dius that determines the 
spatial area over which they respond to scent. The model of Giuggioli et al. ( 201 Ibh 
implicitly identified this perceptive radius with the lattice spacing a. However, this 
limits the model’s flexibility: if the lattice spacing is changed then the model assump¬ 
tions about the animal’s perceptive radius are also changed. Therefore, to ensure our 
model is not constrained by the choice of a, we define the interaction term, C/ {n,m), 
to be a Bernoulli random variable taking value 1 with probability 


1 




(3) 


where j € {1,2}, h{a) is defined so that ah{a) is the perceptive radius of the animal, 
and cU n,m) takes value 0 with probability I — (p. The model from Giuggioli et al.\ 
( 201 Ibl) implicitly had /;(«) = 1. In general, to change the lattice spacing whilst keep¬ 
ing the perceptive radius 5 constant requires setting h{a) = 5 / a, which holds as long 
as h{a) is an integer (see Figure[TJ)). 

The probability ff'inln') of agent i moving to n at timestep m, given that is was 
previously at position n' , is a combination of ^j (n|n') and C/ {n,m), written as follows 


0/{n|n')c/ (n,m) 




dk(n-n') 


if C/(n' -I- 1, m) + Cj (n' — 1, m) 7 ^ 0, 
otherwise, 


( 4 ) 
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where 5^ is the Kronecker delta. 

Equation Q allows us to describe the stochastic algorithm. This is a one-step 
Markov process, so can be fully described by determining the possible states of the 
system at timestep m+l, given the state at time m. Suppose that, for some m, we 
know Ej {n,m) for every n. Suppose further that animal i is at position «,■ at timestep 
m. Then the algorithm is as follows 

1. Calculate £/(n,m + 1) for each n. 

2. Define a categorical distribution taking one of three values n,- — 1, -I- 1 with 

probabilities given by l|n,), ff‘{ni\ni) and + l|n,) respectively. 

These values are the possible future positions of animal i. 

3. Draw a random variable from this categorical distribution and move the animal to 
the position just drawn. 

4. Repeat steps 2 and 3 for each animal in turn. 


2.3 Probability distribution of an agent in a given scent distribution 

To construct a probabilistic master equation describing the above stochastic process, 
we first assume that the evolution of the scent marks can be decoupled from the 
movement of the agent. In other words, we calculate the equation governing a single 
step of each agent’s movement that is true for any fixed, arbitrary scent distribution of 
the other agent. This is a so-called mean-field approximation, that assumes covariates 
between the agent and conspecific scent are small enough to ignore. 

Let U {n,m) (resp. Vin,m)) be the probability of agent 1 (resp. 2) being at position 
n at timestep m and P{n,m) (resp. Q{n,m)) the probability of there being scent present 
of agent 1 (resp. 2) at position n at timestep m. By analysing the probability of moving 
to site n from either site n — I, n, or n + I in one timestep, we eventually arrive at the 
following discrete space-time master equations 


U{n,m+1) =[1 — Q{n + i,m)]^^-U {n — + Q{n + i — 2,m)]+ 

-U{n+ l,m)[l + Q{n + i + 2,m)]^+ 

U{n,m)Q{n + i— l,m)Q{n + i+ l,m), (5) 

V{n,m+ 1) =[l — P{n + i,m)]^^—V{n — l,m)[l -|-P(n-f/ — 2,m)]-|- 

-y(n-|- l,m)[l -|-P(n-|-/-|-2,m)]|-|- 
V{n,m)P{n + i— l,m)P{n + i+ l,m), 


( 6 ) 
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where the following implicit summation notation ( Einsteinl 1916 ) is used 



h{a) — \ 



h{a)— 1 


(7) 


and a is the lattice spacing and the product ah{a) is the perceptive radius of the agent. 

To give some intuition behind equations Q and ([6ll, we focus on equation Q, 
and note that all of the comments in this paragraph hold equally well for equation 
The initial [1 — Q{n + /,ot)] factor in equation (|5]) ensures that there is a low 
probability of moving to position n if there is a high probability of active conspecific 
scent being present at or around position n. The factor [1 + Q{n + i — 2,m)] (resp. 
[1 + Q{n + i + 2,m)]) means that if scent is likely to be present at or around n — 2 
{n + 2) and the animal is at n — 1 (n + 1) at time m then it will be likely to move 
to n at time m + 1. The final summand U{n,m)Q{n + i— 1, m)Q{n + /-f 1,m) means 
that if the presence of scent is highly probable both to the left and right of an animal 
at time m, then it is likely to stay where it is. Notice that if L„t/(n,m) = 1 then 
t/(n, m + 1) = 1 so that probabilities are conserved. 

Let T be the waiting-time between successive jumps. Then equations (|5]l and (|6]l 
rearrange to give 


U{n,m-\- 1) — U {n,m) 



X 




X 


( 9 ) 
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Equation (l8]l can be re-written as follows 

TT^„ I I'l TTf^ 


Ji 

U{n + l,t) — U{n,t) 

1 

1 


a 

a 


( 10 ) 


2a 


1 

2a 

2a 

1 U{n,t)Q{n -f i + l,t)Q[n + i— l,t) — U{n— l,t)Q{n -|- i,t)Q{n + i — 2,f) 
a a 

U{n+ I ,t)Q{n + i + 2,t)Q{n + i,t) — U {n,t)Q{n + i + l,t)Q{n + i— l,t) 


and similarly for equation Taking the limit as a, T —^ 0 and n,m,h{a) ^ °o such 
that D = a^l{2T), x = na, t = mz, ah[a) = 5 in the limit, and writing u{x,t) (resp. 
v{x,t)) for the probability density functions of agent I’s (resp. 2’s) position and p{x,t) 
(resp. q{x,t)) for the probability that agent I’s (resp. 2’s) active scent is present at po¬ 
sition X at time f, we arrive at the following PDE (see Appendix A for a full derivation) 


du 

Tt 


= Z,^[(l-a.l+4D^ 



( 11 ) 


The equation governing the evolution of v{x^t) over time is analogous 


dv 

dt 






dx^ 


dx 



( 12 ) 


Here, p(x,t) and q(x,t) are the locally averaged scent of agents 1 and 2, respectively 

1 

J ^p{x + z,t)dz, (13) 

1 

qix,t) — J^q{x + z,t)dz. (14) 


2.4 Evolution of the scent distribution 

Recall that we gave two different formulae for the scent decay process, equations O 
and (|2]l. Eor the purposes of our mean-field analysis, it is convenient to use equation 
(|2]i. In other words, the probability of scent being present at lattice site n decays by a 
factor of 1 — T/r each timestep of length T. Additionally, when a site is visited by the 
animal, the probability that there is active scent present jumps to 1 with probability 
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The following master equation follows directly from taking the expectation of 
either side of equation (|2|i 

P{n,m+ 1) = (1 — llx)U{n,m)K{a) + (1 — px) [1 —U{n,m)K{a)\P{n,m). (15) 

The probability density version of equation (fTSl) is the limit as a,K{a),x 0 and 
n,m oo of 

p{na,mx + x) =(1 — jlx)u{na,mx)aK{a) + (1 — jix) [1 — u{na,mx)aK{a)]p{na,mx), 

(16) 

such that X = na, t — mx and X — aK{a) /t in this limit. 

Subtracting p{na,mx) from both sides of equation (fThl l. dividing by x and tak¬ 
ing this limit leads to the following ordinary differential equation (ODE) governing 
p{x,t) 


^ p)u-pp. (17) 

We can interpret X as representing the amount of scent deposited over a unit of space 
in a single unit of time. The derivation for q{x,t) is similar and gives 


— =X{l-q)v-pq. (18) 

Analyzing the system of equations (dUl, (HU, ([I7]l and (fTSl) requires choosing 
an appropriate domain and boundary conditions. A simple and biologically realistic 
choice is to assume that agents are confined in a domain [0,L] with zero flux bound¬ 
ary conditions. The boundary conditions could either come about by being confined 
in a valley or on a small island. Alternatively, the conditions could model a situation 
where the rate of migration of animals into the domain is equal to the rate of move¬ 
ment outwards. In other words, the population is assumed to be exhibiting a certain 
spatial and temporal stability. These boundary conditions are given as follows 




dq 

dx 


x=0 




dq 

dx 


= 0, (19) 


x=L 


dx 


[(l-p>]+4 


dx 


x=0 


dx 


[(l-p>]+4 


dx^ 


= 0 . ( 20 ) 


x=L 


The existence of the boundary requires that we need to redefine p{x,t) and q{x,t) 
in the cases where x < 5 and .* > L — 5, as follows 


(jhI^xPix + z,t)dz ifx<d, 

p{x,t)= < ^!^sp{x + zd)dz if5<x<L-5, 

+ ifx>L-5, 


( 21 ) 
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q{x,t) 


h!^5^{^ + Z,t)Az 


ifx < 5, 

if 5 < ;c < L — 5, 
'\f x> L — 5. 


( 22 ) 


In addition to the boundary conditions, it is necessary to impose integral condi¬ 
tions on the initial probability distributions u{x,Q) and v(x,0), to ensure that proba¬ 
bility is conserved. In other words 

f u{x,Q>)&x= [ v(.x,0)cbi: = 1. (23) 

Jo Jo 

A consequence of equations ( fT9b and (l20l) is that the time-derivative of Jq u{x^t)Ax is 
zero. Therefore the initial conditions from equation (l23t imply that probabilities are 
conserved at every point in time, i.e. 


/ u{x^t)dx = / v(.r,f)dx=l. 

Jo Jo 


(24) 


2.5 A dimensionless version of the model 

To minimize the number of model parameters, we re-write equations (fTTI) . (fT^ . (fTTI) . 
and (fTsT l. using the following dimensionless parameters 

X ^ tD jxL D 

u = Lu,v = Lv,x=-,t = -px,m =e = (25) 
L A LA 

Dropping the tildes over the letters to ease notation, we arrive at the following di¬ 
mensionless system of equations, which will be the object of study for the rest of this 
paper 


du 

dt 


dx^ 


[n-m+iA 



dv 

dt 


dx^ 


[(l-/)y]+4 


dx 



e 


dt 


(1 — p)u — mp, 


(26) 


(27) 


(28) 


e 


dq 


(1 — q)v — mq. 


(29) 
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3 Territorial patterns 


We define a territorial pattern to be a non-trivial steady-state solution to equations 
(l26l l-(l29ll. These are found by setting to zero the left-hand sides of equations (l26l l- 
( |29] |. Setting equation (l28l l (resp. equationl29ll to zero enables the steady state solution 
of p{x,t) (resp. q{x,t)), denoted by p*{x) (resp. q*{x)), to be written in terms of the 
steady state solution of u{x,t) (resp. v{x,t)), denoted by u*{x) (resp. v*(jic)) as follows 


p*{x) 


u*{x) 

m + u* (x) ’ 


(30) 


q*{x) = 


\x) 


m + v*{x)' 


To ease notation, we will henceforth drop the asterisks. By setting equations 
(l27l l to zero and integrating with respect to x, we have that 


dr 


{{1 - q[v{-),x]^)u{x)} +4 


dq 

dr 


ix) 


= CI, 


(31) 

and 

(32) 


dr 


{{I -p[u{-),xf)vix)} +4 


dr 


v(x) 


= C2, 


(33) 


for constants ci and C 2 - The boundary conditions given by equations ( fT^ and (l20l) 
imply that ci = C 2 = 0. 

We use the notation p[u{-),x\ and ^[v(-),r:] to emphasize the fact that p and q are 
functionals. That is, they map the functions m(-) and v(-), respectively, to the interval 
[0,1]. These functionals are given by the following formulae 


p[u{-),x] = 




1 ujx+z) 

1 rS ujx+z) 

^ 5 m+ujx+z) ^ 


ifx<d, 

if 5 < r < 1 — 5, 

1 rl-x ujx+z) 

: J — 5 m+}i(r+7)^ It X > i O, 


l—x+d-> — 5 m+ujx+z)^ 


-L. f _ 

x+5 u—x m+vjx+z) 


vjx+z) 


dz, 


vjx+z) 


■dz 


1 r _ 

2^ t —5 m+vjx+z^ 

1—r+t) f—S m+vjx+z) 


ifx<d, 

if 5 < r < 1 — 5, 


(34) 


(35) 


1 ifx>l-5, 


In sum, as well as equations (l34l) and dTSl) . we have the following system of 
equations, whose non-constant solutions correspond to territorial patterns 


dr 




dq 

dr 


(r) 


= 0 , 


(36) 


dr 


{{I - p[ui-),x]^)v{x)} +4 


dr 


v(r) 


= 0 , 


(37) 
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(38) 



(39) 


3.1 Territorial patterns with only local interactions 

We first examine the case where 5 —?> 0 so that agents only respond to scent at the 
exact position where they are situated. This means equations (l36l l and dJTl i become 


^{{\-q{xf)u{x)}+A- ^u{x) =0, 


(40) 


i{(i-,(xm+4fe,w =0. 


(41) 


The limit 5—^0 means that the functionals p[u{-),x\ and have been replaced 

by functions p(x) and q{x), which makes analysis tractable. To ease notation, we 
hencefore drop the explicit dependence of the functions m, v, p, and q on x. 

By substituting equations dMt and (l39]) into (l40l i and (HTT) . the following system 
of ODEs for the steady state solution of (m, v) is found 


Au = 0, 



(42) 


The system of ODEs in equation (l42l l is simple enough to analyze mathematically. 
The results of this analysis are summarized in the following 

Theorem 1 7. No scent decay. If m = 0 then p(x) = q{x) = 1 and u(x), v(x) can 
take any value. 

2. Positive scent decay. Ifm>0 then there are no non-constant solutions to equa¬ 
tion (ES- Hence no territorial patterns can form in this case. 


Proof See appendix B. 
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3.2 Territorial patterns with non-local interactions 


In the case where 5 > 0, equations (l36l) -(l39b give a system of integral-ODEs, so are 
harder to analyse a nalytically. Instead, we solve them numerically using the method 
of false transients dMallinson & de Vahl Davis , 1973h . This involves solving equa¬ 
tions (I26l l- (l29l l forward in time until the solution is unchanging. 

Our algorithm uses a forward-difference approximation for time and a central 
difference approximation for space. We divide the interval [0,1] into 1,000 equal, 
non-intersecting, sub-intervals of length 0.001. We iterate finite-difference versions 
of equations (|2^ -(|29]| using timesteps of 0.01, until all of the u{x,t) or v{x,t) values 
in all of the sub-sections are increasing by less than 10^^ over each timestep. The 
initial conditions have all of u{x) concentrated on the sub-interval [0.25,0.251) and 
all of v(.r) on the sub-interval [0.75,0.751). This means u{x) and v(.r) are zero outside 
the sub-intervals [0.25,0.251) and [0.75,0.751) respectively, and each integrate to 1 
over [0,1]. 

Numerical analysis shows that patterns emerge from this system corresponding 
to two territories: u{x) on the left and v(x) on the right (figure |2^,b). Notice that a 
larger scent averaging radius leads to wider overlap of the probability distributions, 
meaning that the perceptive scale of the animal plays a large role in the territorial 
patterns that emerge. 

These can be compar ed with the terri t ories th at form in the original IBM with 
the interaction rules from Giuggioli et al\ ( 201 lb ). Although there is some qualita¬ 
tive agreement, the patterns generated by the IBM are still quite different to the PDE. 
In the IBM, at any point in time, there is a border between the two territories. This 
border fluctuates about the central point, typically much slower than the movement 
of the agent. Each agent is free to move within its territory borders. Consequently, 
the probability density of both agents combined (m + v) ends up being roughly uni¬ 
form (figure|2};). This does not happen in the mean field approximation studied here. 
Indeed, the value of m -f v appears to be lower in the middle of the terrain. Since this 
is just an artifact of the assumptions made in using the PDE limit, it is necessary to 
be cautious when inferring biological lessons from such pattern features. 


4 Investigating pattern formation via linear analysis 


A common technique for examining whether patterns spontaneously form in a dy¬ 
namical system is to linearize the syst em abou t the un iform steady state and examine 
the resulting dispersion relation, e.g. Murray ( 2002h chapter 2. Eor our system, the 
uniform steady state is 


{us,Vs,Ps,qs) 



1 -f m ’ I + m J 


(43) 


That Mi = Vj = 1 arises from the integral conditions (equation |24|. The values for ps 
and qs then follow from equations dMt and (l39l l. 
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Fig. 2 Numerical steady state solutions of the model. Solid red (resp. black) lines denote values of z<(x) 
(resp. v(a:)), whereas dotted red (resp. black) lines show values of p{x) (resp. q[x)). In both panels, m — 0.4 
and e = 0.01. In panel (a), 5 = 0.1, whilst panel (b) has 5 = 0.01. Notice that a larger the scent-averaging 
radius, 5, gives a larger overlap between temtories. Panel (c) compares steady states of simulations of the 
original IBM (dashed lines) with numerical results from the PDE approximation (solid lines). As in panels 
(a) and (b), dotted red (resp. black) lines show values of p{x) (resp. q(x)). Here, 8 — 0.01, m — 0.4, and 
£ — 0.01. This corresponds, in the IBM, to Tas — 500 and N — 100. 




Fig. 3 Dispersion relations. Panel (a) show the dispersion relation for the dynamical system in equations 
126I29L where m = 0.1 (solid line), m — 0.5 (dotted), m— I (dashed), and m — 3 (dot-dashed). We set 
£ — 0.01 and 5 = 0.01 throughout. Panel (b) shows the same dispersion relations, but this time the animals 
respond only to the scent density at the particular point in which they reside, i.e. we take the limit 5 —> 0. 
This system is given in equations 1281 . 1291 . j461 and 1471 . The values of £ and m are identical to those in 
panel (a). 


Letting w = {u,v,p,q) = {u — Us, 
to give the linearized system 


du 

dt 

dv 

dt 


1 - 

1 - 


1 


(1 +m) 
1 

(1 +m) 


dt 

dq 


m 

-tM 

e(l +m) 

m 

e(l +m) 


V — Vs,p — ps,q — qs), use equations (I26H29I ) 


d'^ii 

9 

1 


1 +m 

dx^^ 

9 

1 


1 +m 


1 +m ^ 


1 + m ^ 


d^q 
dx^ ’ 
d^p 
dy?- ' 


( 44 ) 
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Searching for solutions of the form w = (mo, vo,po,?o)exp((7f + ikx), we obtain the 
following eigenvector equation 


Aw = (7W 


/ 


1 

(T+mp' 





A = 


V 


0 

m 

e(l+m) 

0 


0 


1 

(T+mp 



0 




m 

£(l+m) 


0 

-2[2-t^] I sin 5^ 

_ 1+m 

e 

0 


-2[2-T^]|sin5^^ 

0 

0 

_i^ / 

(45) 


The dispersion relation is given by plotting the real values of <7 as a function of the 
wave number k, wherever det(A — al) — 0. As shown in Figure[3^, patterns can form 
for a finite range of wavelengths as long as m < 1; that is, as long as the scent decay 
is not too rapid. 

We can gain biological insight by relating this result back to the underlying IBM. 
Recall that m = jiL/X (eauationlZSll. Recall also that A is the limit of aK{a)/z. In the 
original lattice model, where K{a) = 1, aK{a)/z is simply the speed of the animal. 
Then m < 1 if and only if the time it would take a freely moving animal on the lattice 
to traverse the whole terrain is less than the characteristic timescale for scent-mark 
decay 1/^. 

The dispersion relation changes somewhat if we examine the case where 5—^0, 
so that animals only respond to scent in the exact place that they are located at any 
point in time. In this case, p and q are replaced by p and q respectively, so that 
equations (|26]) and dZTl i are replaced by 


du 

dt 


dx^ 


[d-,=)„]+4|- 



(46) 


dv 

dt 


dx^ 





(47) 


Equation (l45l l becomes 


/ 


A = 


(T+mp 


-1 


0 

m 

£(l+m) 

0 


0 


(l+m)‘ 


■- 1 


k ^ -2 [2 - 


e(l+m) 


^1 
1+m J 

1+m 

£ 

0 


-2[2-Ti;]^"\ 
0 
0 

_ 1+m 

£ / 


(48) 


which is the limit as 5 —?► 0 of equation (l45T l. The corresponding dispersion relation is 
given in Figure [2t. Here, for 0 < m < 1, O’ is an increasing function of k, indicating 
that the steady state is unstable but arbitrarily large wave numbers grow fastest. In 
other words, this is an ill-posed problem. 
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5 Discussion and conclusions 


We have shown how stable territorial patterns can form purely from a conspecific 
avoidance mechanism, without requiring any attractive potential. Our model is con¬ 
structed by taking the continuous space-time limit of a discrete lattice model. There¬ 
fore it can be rigorously linked to the underlying movement and interaction processes. 
We have demonstrated that patterns will only form if the scent marks last for a suffi¬ 
ciently long time. If they decay too quickly, i.e. m > 1, the territorial structure breaks 
down. This can be interpreted as saying territories can only emerge if the animal is 
able to patrol its territory faster than the scent marks decay. 

Similarly, patterns will only form reliably if the animals react to the averaged 
scent density across the local vicinity of the animal. From a biological perspective, 
an animal will always have a perceptive radius over which it will react to scent. 
Therefore this spatial averaging is implicit in the system being modeled. As such, 
our study demonstrates the importance of ensuring that the mathematical limiting 
process, moving from discrete to continuous space, does not destroy a key feature of 
the underlying biology. Our procedure for performing this limiting process has the 
potential for broad application, since there are many examples where the discreteness 
of ecological interaction s is known to be an important feature of the modeling process 
dPurrett & Levinlll994 . 

The model is deriv ed from an individual-based model, previ ously studied using 


stochasti c simulations dGiuggioli et oDi 201Jj Pottsef^ 


_ _ _ _ l2Q12h . As noted in recent 

reviews ( Giuggioli & Kenkm 2014 : Potts & Lewis . 2014 ), one of the advantages of 
this approach is that it gives a clear delineation between the related notions of ‘home 
range’ and ‘territory’. The territor y of an anim al is defined as the area containing 
active scent marks of the animal dSurl . 1943h . Therefore, in the model presented 
here, p{x,t) and q{x,t) can be considered the probabilities of position x being part of 
the animals’ territories at time t. 

_ O n the other hand, the home range of an animal is its utilization distribution 

dBurlill943h . Therefore u{x,t) and v(x,t) can be considered as the home ranges of the 
animals at time t. The utilization distribution of an animal is typically m uch easier to 


measure in the field than the fluctuating locations of the territory border (IPotts et al. 


20121) . In our approach, the concepts of territory and home range are related by rather 
simple steady-state equations (l30l l and (OTT i. This gives an explicit way to calculate 
the probable location of a territory border, given data on its home range. 

A key reason for studying PDF limits of IBMs is to provide mathematical analysis 
of the conditions under which patterns may form, rather than relying on empirical ev¬ 
idence from computer simulations. However, as shown here, patterns that form from 
numerically solving the PDF may to be quantitatively different from those formed 
by simulating the IBM. Therefore, if such PDF models were fitted to data on real 
systems, it is important for the user to check that the PDF results are not signifi¬ 
cantly different to those given by the IBM. Otherwise, there is a danger of making 
incorrect inferences about biological patterns, that may merely arise as artifacts of 
the mean-field approximation and/or limiting procedure. 

Models such as ours could be of use in analyzing territory formation when there 
is no reason to believe the animals have any fidelity towards particular locations, or 
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where these locations are not known, e.g. iBateman et al\ (12015h . Th ough memory 


processes have recentl y been invoked to explain pattern formation (iBriscoe et al. 


2002 ; Moorcrofl 2012 ), it is unclear how to find out what is going on inside the minds 


of the animals using current science. This makes conjectures about memory difficult 
to falsify. C onspecific avo i dance mechanisms, on the other hand, can be measured 
directly, e.g. Arnold et al\ ( 2011 ). Therefore our model of territorial emergence has 
the potential to be parametrized from empirically measured interaction mechanisms. 
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Appendix A 


Let u{x,t), v{x,t), p{x,t) wdq{x,t) be the density functions corresponding to U {n,m), 
V{n,m), P{n,m) and Q{n,m) respectively, where x = an and t = niT. First note the 
following limit as a ^ 0, k{a) ^ ak{a) ^ 5 


1 

2k{a) - 1 


k{a)— 1 

L Q(n + i,m) 

i=\—k{a) 


1 

2k{a) - 1 


*(")-! J nS 

Y q{x + ai,t)a^ — / q{x + z,t)dz. 
i=l k{a) 25 J-S 


(49) 


Using the definition of q{x,t) given in equation (O, and writing equation (fTOl l down 
in terms of the density functions, we have 


i{x,t + t) — u{x,t) 



u(x + a,t) — u(x,t) 

u{x,t) — u{x — a,t) 

\ a 

a 

a 


1 

2a 
1 ^ 
a 


A ( , .q{x + 2a,t)-q{x,t) q{x,t) - q{x-2a,t) 

4u(x + a,t) --- 4u(x— a,t) - 


2a 2a 

u{x,t)q{x + a,t)q{x — a,t) — u{x — a,t)q{x,t)q(x — 2a, t) 


u{x + a,t)q{x + 2a,t)q(x,t) — u{x,t)q{x + a,t)q{x — a,t) 


(50) 


We keep x constant in the limit as a —0,n —oo. Taylor expanding the right-hand 
side about x, assuming a is arbitrarily small, gives the following expression 


u{x,t 4-x) — u{x,t) a^ f d^u 


2x 1 dx^ 


du dq d^q 

4 -^+4m—^ 

dx dx dx^ 




dq dq d^u 


- 2u --—I- ^ q^ 4~ ^ -- q 

dx dx dx^ dx dx ^ 


du dq _ 


+ 0{a) 


(51) 


In the limit as a, T —?► 0 such that a^ j (2t) — D, this simplifies to give equation (fTTIi . 
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Appendix B 

We look for solutions to equation (l42l) in two cases: det(A) ^ 0 and det(A) = 0. For 
det(A) 7 ^ 0, one solution is to have u = 0, implying that u{x) and v{x) are constant 
functions so territorial patterns do not form. 

Otherwise, suppose that det(A) ^ 0, du/dx = 0 and dv/dx ^ 0. Then the following 
equations hold 


2u{2m + v) = 0, (52) 

{m + 2u)(m +u) = 0. (53) 

Equation (l52l i implies m = 0 or v = —2m. However, if m = 0 then equation (l53T l would 
imply m = 0, which contradicts det(A) ^ 0. Furthermore, if v = —2m then v < 0, 
which contradicts the fact that v(x) is a probability density function. In conclusion, 
if det(A) 7 ^ 0, we cannot have du/dx = 0 and dv/dx 7 ^ 0. Similarly, if det(A) 7 ^ 0, 
we cannot have dv/dx — 0 and du/dx^O. Therefore the only possible way for non¬ 
constant steady states to arise is if det(A) = 0. 

Lemma 1 Ifdet{A) = 0 then there are two possibilities. 

1. No scent decay. Ifm = 0 then p(x) = q(x) = 1 and u{x), v{x) can take any value. 

2. Positive scent decay. If m > 0 then, for each x G [0,1], there are finitely many 
possible values for u{x) and v{x), one of which is u(x) = v{x) = m. Furthermore, 
all solutions other than u(x) = v{x) = m have u{x),v{x) f m. 

Proof If m = 0 then det(A) = 0. Furthermore, by equations (l30l l and OTT i. we have 
p{x) = q{x) = 1. Therefore dp/dx = dq/dx = 0 so that equations (l32T i and (l3Jt hold 
regardless of the values of u{x) and v(x), proving part 1 of the lemma. Indeed, in the 
time-dependent PDFs (I28I29I46I471 |. if an initial condition of p(x,0) = q{x,0) = 1 is 
given and m — 0 then dp/dt — dq/dt — du/dt — dv/dt = 0 so uixf) = u{x,Q) and 
v{xf) = v(x, 0 ) remain unchanged for all times t. 

Now suppose mfO. For notational ease, we drop the explicit dependencies of u 
and V on X for the rest of this proof, noting that they always refer to the steady states. 
Then the equation det(A) = 0 implies the following polynomial holds 

(m -f 2v) (m + 2u) {m -f u) {m -f v) = 4Mv(2m + u) (2m -f v). (54) 

Equation (l54l i can be rearranged to give 

u^{2m — 2v) -f u{3m^ — 7mv — 2v^) + (m^ -f 3m^v-|-2mv^) = 0. (55) 

Clearly m = u = v satisfies equation (l55l l. Furthermore, it follows from equation (l55l l 
that u — m if and only if v = m. Hence any solution other than u = v = m has both 
u f m and v f m. 

Differentiating equation (l54l i with respect to x, we find 

du dv7mu — 3m^ — 4mv + 2u^ + 4uv 
dx dx 3m2 — 7mv-|-4mM —2v^ —4 mv 


( 56 ) 
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Using equation (l5^ and the top line of the vector equation (l42li . together with our 
assumption that du/dx,dv/dx ^ 0, we find that 


(7mM — — Amv + 2u^ + 4uv) (m + 2v) (m + v) + 

[Amu + 2mv) (3m^ — lmv + Amu — 2v^ — Auv) = 0. 


(57) 


Proving part 2 of the lemma requires applying Bezout’s Theorem ( Fultoni 1969ll 
to equations (l55l l and (fSTl l. Bezout’s Theorem states that if two projective plane curves 
are zeros of polynomials with no non-constant greatest common divisor, then the 
curves intersect at finitely many points. The polynomials on the left-hand sides of 
equations (l55l l and (fSTl i are homogeneous in three unknowns, therefore equations (fSSl) 
and (fSTl i describe curves in the real projective plane. Thus, to prove part 2 of Lemma 
[T] it suffices to show that these two polynomials have no non-constant common factor. 

Let /(m, M, v) be the polynomial on the left-hand side of equation (l55l l. Since this 
is quadratic in u, it written as precisely one of the following two possible decompo¬ 
sitions; 


f{m,u,v) = [ai{m,v)u + bi{m,v)][a 2 {m,v)u + b 2 {m,v)]c{m,v), (58) 


or 


f{m,u,v) = a{m,u,v)l5{m,v) 


(59) 


where ai (m, v), 02 ( 111 , v), b\(m,v), b 2 (m,v), and j3 {m,v) are polynomials, and a{m,u,v 
is an irreducible polynomial. By solving equation (fSST l in terms of u, we find that 


7mv + 2v^ — 3m^ ± \/A(v — J+m) (v — 7 _m) (v — t?+m) (v — j?_m) 

A(m — v) 


(60) 


where 



-ll-4v^±2j^-f30V3 


t?± = -y + y3±i7l5(15-8V3). 


(61) 


Therefore the numerator of equation dhOl l is not a polynomial, so the decomposition 
given in equation (fSHT l cannot hold. 

It follows that f{m,u,v) = a{m,u,v)P{m,v) where a{m,u,v) is irreducible and 
I5{m,v) is the greatest common divisor of the coefficients of m" in f(m,u,v) for 
n = 0,1,2. These coefficients are 2m — 2v, 3m^ — Imv — 2v^ and m^ -f 3m^v + 2mv^ 
(eauationl55]). Since m — v does not divide 3m^ — Imv — 2v^ or m^ + 3m^v A- 2mv^, it 
follows that p (m, v) is a constant. Hence f(m,u,v) is irreducible. 

Since f{m, u, v) is both irreducible and not a constant multiple of the polynomial 
in equation (fSTl i. there is no non-constant greatest common divisor of the polynomials 
in equations (l55l l and (fSTT i. The proof of part 2 then follows from Bezout’s Theorem. 

□ 
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Fig. 4 Possible values of u,n — u(x)/m and — v(x)/m. The black curve denotes solutions to equation 
<62), whereas the grey curve shows solutions to equation ii). There appears to be only one crossing-point 
for positive real values of both and v^, which is where = v„, = 1, so that u — v — m. 


Lemma[T]enables us to prove Theorem[T]from Section[3] as follows. 

Proof of Theorem\J] Part 1 of Theorem [T] is identical to part 1 of Lemma[T] To show 
part 2, note that classical solutions must be continuous. Lemma[T]states that there are 
only finitely many possible values of u and v. Therefore any classical solution must 
be constant. □ 


Note 1 Numerical analysis suggests that u = v = mis the only positive real solution 
(Fig.|4|i. Since we are interested in the case mfQ,we set Vm = vjm, = ujm and 
assume v f^m. Then equation (l55l l implies 


_ 7Vm + 2v^ - 3 ± •\/(3 - Ivm - 2v2 )2 _ 8(1 - Vm)(l + 3Vm + 2v2 ) 

Urn — \ 

4(1 - Vm) 

Furthermore, equation (l57l i rearranges to give as another two-valued function of 


4v)^ + 4v^ + 3vm+ 19 

Mm =--± 


8v2-|-4vm-36 

\/(4Vm T4vm -f 3vm + 19)^ ~ 4(4v^ + 2vm — 18)(8v^, + 18v^ + 13vm + 3) 

8v2+4vm-36 


( 63 ) 
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The black curve in Fig. |4]has an asymptote at Vm = 1, where the denominator of 
the right-hand side of equation (l62]) tends to 0. The grey curve has an asymptote at 
Vm = (v/^— l)/4, where the denominator of the right-hand side of equation (l6^ 
tends to 0, so the two curves do not cross at values of higher than those shown in 

Fig. a 
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